Behavioral Segmentation

Behavioral analysis has become one of the most important advancements in the movement ecology literature. Due to the importance of the internal state in individual movement decisions, analyses such as these offer an opportunity to understand the motivations underlying space-use more clearly than most of the broader scale analyses we’ve seen thus far. We are going to begin by a simple example of First Passage Time (FPT), then consider more complex methods such as Hidden Markov Models (HMM) and Behavioral Change Point Analysis (BCPA).

library(tidyverse)
library(adehabitatLT)
library(sf)
library(moveHMM)
library(bcpa)

#load albatross data
albatross <- read_tsv("data_files/Study2911040", comment="##", col_types = "Tnncc") %>%
  as.data.frame() %>%
  na.omit()

#add UTM coordinates
coords <- sp::SpatialPoints(data.frame(albatross$location_long, albatross$location_lat), proj4string=CRS("+proj=longlat +ellps=WGS84")) 
UTM_coords <- spTransform(coords, CRS("+proj=utm +south +zone=16 +ellps=WGS84")) # reproject to UTM (meters)
albatross$utm_long <- UTM_coords@coords[,1]
albatross$utm_lat <- UTM_coords@coords[,2]

First Passage Time (FPT)

First passage time (FPT) measures the length of time an animal spends within an circle of a given radius and is commonly used as a proxy for area-restricted search behavior (Fauchald & Tveraa 2003). We first calculate the FPTs for a range of radii, and plot the variance in FPT vs. radius to identify a ‘characteristic scale’ of restricted search behaviour where the variance is maximized. FPT is easily calculate in the adehabitatLT package, so first we need to convert our data into an ltraj object. Let’s start with an FPT analysis on our Galapagos Albratoss tracking data.

# create 'ltraj' object
alba_ltraj  <- as.ltraj(xy=cbind(albatross$utm_long, albatross$utm_lat),
                        date=albatross$timestamp,
                        id=as.factor(albatross$individual_id),
                        proj4string = CRS("+proj=utm +south +zone=16 +ellps=WGS84"))

# Calculate fpt in hours for one individual (#2911074) in radii of 0-400 km in 10-km increments
albatross_fpt <- fpt(alba_ltraj[10], radii = seq(0, 400000, 10000), units="hours")

# Calculate the variance for log(FPT) for radii and find the peak radius
fpt_var<-varlogfpt(albatross_fpt) 

There appears to be a peak at very small spatial scales (1 km) and a secondary peak at large spatial scales (about 120 km). Let’s investigate the behavioral patterns at large spatial scales.

# plot fpt over time
plot(albatross_fpt, scale=120000, xlab=c("date")) #scale = radius to plot at

# add fpt at characteristic scale to tracking data
albatross10 <- albatross %>% filter(individual_id=="2911074")
albatross10$fpt <- albatross_fpt[[1]]$r12

# plot track by fpt
ggplot(albatross10) +
  geom_point(aes(x=location_long, y=location_lat, color=fpt)) + 
  scale_color_gradientn(colors=rev(topo.colors(10))) +
  theme_bw()

Here we can pick out based on FPT high-use areas around the albatross’s breeding area on the Galapagos Islands to the East, and foraging areas during the non-breeding season off of mainland South America to the East.

Hidden Markov Models (HMM)

Let’s now fit a Hidden Markov Model to the same albatross track. Note that we can fit these models at once for all individuals in the population, but to make our demonstration run faster we will just use the one individual. The first step in the model building procedure is using the moveHMM::prepData command. We will specify our data, the type of projection (in our case ‘LL’ for latlong coordinates), and the names of the columns representing our coordinates (location_long and location_lat):

data <- prepData(albatross10, type="LL", coordNames=c("location_long","location_lat"))

## Warning in prepData(albatross10, type = "LL", coordNames =
## c("location_long", : There are 168 missing covariate values. Each will be
## replaced by the closest available value.

head(data)

##        ID        step     angle         x         y           timestamp
## 1 Animal1 0.005007388        NA -89.62092 -1.389501 2008-06-23 17:58:44
## 2 Animal1 0.011744984  2.968880 -89.62097 -1.389518 2008-06-23 18:01:59
## 3 Animal1 0.013633401  3.124203 -89.62086 -1.389495 2008-06-23 19:30:39
## 4 Animal1 0.007631059  2.749848 -89.62098 -1.389520 2008-06-23 21:00:49
## 5 Animal1 0.010553185  3.045311 -89.62092 -1.389533 2008-06-23 22:30:24
## 6 Animal1 0.006036250 -3.121484 -89.62101 -1.389506 2008-06-24 00:00:56
##   individual_id  tag_id utm_long utm_lat      fpt
## 1       2911074 2911108 208339.6 9846256 8.231712
## 2       2911074 2911108 208335.0 9846254 8.231712
## 3       2911074 2911108 208346.4 9846257 8.231712
## 4       2911074 2911108 208333.1 9846254 8.231712
## 5       2911074 2911108 208340.6 9846253 8.231712
## 6       2911074 2911108 208330.4 9846256 8.231712

Now we have an object (data) with three new variables: ‘ID’, ‘step’ (step length in coordinate units) and ‘angle’ (in radians; i.e., ranging from -pi to pi). We can also visualize the path and step length and angle distributions using the plot command. This will also give us time series plots of the step length and turning angle.

plot(data)

Now it is time to fit an HMM to the data. To do this, we will use the moveHMM::fitHMM command. This is a pretty complex function, however, that requires quite a few inputs to make it run smoothly. Ultimately, our goal is to build a two-state model that relates the behavioral state to the turn angle and step lengh distributions.

Unless you have ancillary behavioral information (e.g. dive depth, accelerometry, etc.) it is best practice to start with nbStates=2 which specifies that we want to fit a two-state model to the data. Then we will offer a formula to guide the calculation of the function regarding transitions between states. A formula=~1 indicates no covariate effect. If we had a covariate, like chlorophyll or NDVI, we would use formula=~chlorophyll. Next we want to define the distrbutions that we want to use to characterize both the step lengths and turning angles. We are going to use a gamma distribution for the former (stepDist="gamma") and a vonMises distribution for the latter (angleDist="vm"). That takes care of that, but we’re still not ready to run our function just yet. We want to define some initial values for the state-dependent functions so that the optimization algorithm has a starting point. In this case, the initial parameters should be specified in two vectors, stepPar0 (for the step distribution) and anglePar0 (for the angle distribution). The exact parameters needed for each can be found in the documentation, but for a gamma distribution, we will need a mean, SD, and zero-mass (the latter only if there a step lengths equal to 0) and for the vonMises, we will need a mean and concentration parameter. This vignette has a nice explanation for choosing starting parameters, but a good place to start is looking at the empirical step length and angle distributions.

mu0 <- c(20,40) # step mean (two parameters: one for each state)
sigma0 <- c(20,5) # step SD (two parameters: one for each state)
(stepPar0 <- c(mu0,sigma0))

## [1] 20 40 20  5

angleMean0 <- c((pi/6),0) # angle mean 
kappa0 <- c(1,1) # angle concentration (two parameters: one for each state)
(anglePar0 <- c(angleMean0,kappa0))

## [1] 0.5235988 0.0000000 1.0000000 1.0000000

Now we can fit our model with the moveHMM::fitHMM command:

m <- fitHMM(data=data, nbStates=2, stepPar0=stepPar0, anglePar0=anglePar0, formula=~1, verbose=T)

## iteration = 0
## Step:
##  [1] 0 0 0 0 0 0 0 0 0 0 0
## Parameter:
##  [1]  2.9957323  3.6888795  2.9957323  1.6094379  0.8660254  1.0000000
##  [7]  0.5000000  0.0000000 -1.5000000 -1.5000000  0.0000000
## Function Value
## [1] 7114.632
## Gradient:
##  [1]  5928.9975847  -193.1542842 -5216.8284156   -42.9859826    37.4064120
##  [6]   -34.9566762   288.7923347    -1.6976564   192.0960003   -13.4982887
## [11]     0.5000002
## 
## iteration = 107
## Parameter:
##  [1]   1.61836319   3.73474390   2.21667993   3.10016224   0.59565035
##  [6]   9.18705903  -0.02746546   0.05141144  -3.40944914  -1.57911558
## [11] -14.03357793
## Function Value
## [1] 5138.606
## Gradient:
##  [1]  1.700564e-03 -4.578226e-04 -1.187807e-03  1.106005e-04 -1.273293e-04
##  [6] -7.632697e-05 -2.346496e-04 -1.291482e-04 -1.493840e-05 -9.388017e-05
## [11]  7.777016e-07
## 
## Relative gradient close to zero.
## Current iterate is probably solution.

Fortunately, that doesn’t take too long even though there are some pretty intense calcualtions going on in the background. This is primarily because we are fitting relatively few data points (only 1326 in total). We can take a look at the resulting model m:

print(m)

## Value of the maximum log-likelihood: -5138.606 
## 
## Step length parameters:
## ----------------------
##       state 1  state 2
## mean 5.044826 41.87730
## sd   9.176813 22.20155
## 
## Turning angle parameters:
## ------------------------
##                  state 1     state 2
## mean          -0.0460774 0.005596013
## concentration  0.5962832 9.187202882
## 
## Regression coeffs for the transition probabilities:
## --------------------------------------------------
##              1 -> 2    2 -> 1
## intercept -3.409449 -1.579116
## 
## Transition probability matrix:
## -----------------------------
##           [,1]       [,2]
## [1,] 0.9679985 0.03200146
## [2,] 0.1709208 0.82907923
## 
## Initial distribution:
## --------------------
## [1] 9.999992e-01 8.040706e-07

This output has all sorts of interesting information for us. The first thing we can see is a log-likelihood value. Good to know, but not especially meaningful by itself. Next, we have the step length parameters. The model has explored parameter space for the mean and SD parameters and returned optimal values of each for both of our behavioral states. We can see right off the bat that the mean step size of State 1 is quite a bit smaller than that of State 2, so we have some idea about what kind of activities may be occuring in each. We may have something like Area Restricted Search (ARS) during State 1 and more directional movement during State 2. The next section defines the turning angle parameter estimates. Next, we can see the regression coefficients for the simple formula we set up to calculate state transition probabilities.

We can also use the moveHMM::plot command to visualize all of these things, beginning with the distributions of step lengths and turning angles in the two states, illustrating the transition proababilities between states, and then showing each of the four paths with each point assigned to the most likely state.

plot(m)

## Decoding states sequence... DONE

We’ve officially built our first hidden Markov model! Let’s see some ways that we can use the model outputs. The first is to ‘decode’ the behavioral states along the paths. This was done for us when we plotted each track above, but if we wanted to see the most likley states for each point, we could use the moveHMM::viterbi command, which uses the Viterbi algorithm to predict the most likely sequence of states that generated these paths. From this, we can determine the proportion of time that the albatross spent in one state versus the other:

states <- viterbi(m)
prop.table(table(states))

## states
##         1         2 
## 0.8484163 0.1515837

It turns out that this animal was in the ARS state (State 1) nearly six times as frequently as they were in the directional movement state.

If we wanted to get a little bit more information of the probabilities of the animal being in a particular state at a given time, we could use the moveHMM::stateProbs command on our model. In this case, rather than extracting one most likely state for each point, the actual probabilities of both states are displayed for each point.

state.probs <- stateProbs(m)
head(state.probs)

##      [,1]         [,2]
## [1,]    1 1.398323e-18
## [2,]    1 1.649479e-19
## [3,]    1 2.332904e-19
## [4,]    1 6.951051e-20
## [5,]    1 1.059858e-19
## [6,]    1 1.640253e-20

We can also visualize these probabiities using moveHMM::plotStates:

plotStates(m)

## Decoding states sequence... DONE
## Computing states probabilities... DONE

Now let’s see how our model performs with an ancillary variable like fpt. Then, we can use the moveHMM::AIC command to compare the likelihood values (or more accurately, the information criterion values derived from the liklihoods) and determine which model performed better.

m2 <- fitHMM(data=data, nbStates=2, stepPar0=stepPar0, anglePar0=anglePar0, formula=~fpt, verbose=T)

## iteration = 0
## Step:
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0
## Parameter:
##  [1]  2.9957323  3.6888795  2.9957323  1.6094379  0.8660254  1.0000000
##  [7]  0.5000000  0.0000000 -1.5000000  0.0000000 -1.5000000  0.0000000
## [13]  0.0000000
## Function Value
## [1] 7114.632
## Gradient:
##  [1]  5928.9975847  -193.1542842 -5216.8284156   -42.9859826    37.4064120
##  [6]   -34.9566762   288.7923347    -1.6976564   192.0960003 26917.9309562
## [11]   -13.4982887 -1058.7997640     0.5000002
## 
## iteration = 157
## Parameter:
##  [1]    1.69260285    3.78060846    2.29142271    3.08342031    0.61985305
##  [6]   11.88148523   -0.02570511    0.02565875   -2.40117704   -0.01479218
## [11]   -2.19361668    0.01857612 -117.79988317
## Function Value
## [1] 5123.122
## Gradient:
##  [1]  2.348154e-04 -8.468006e-05 -2.099668e-04 -2.477689e-05  5.820766e-05
##  [6]  7.654722e-07  4.638423e-05  2.546585e-05  9.848029e-06  1.091394e-04
## [11] -4.560706e-06 -5.638867e-04  0.000000e+00
## 
## Relative gradient close to zero.
## Current iterate is probably solution.

AIC(m, m2)

##   Model      AIC
## 1    m2 10272.24
## 2     m 10299.21

The lower AIC value for m2 indicates that we got additional information that was more beneficial than the additional cost of adding a another parameter. Now we can take a look at our more accurate model:

plot(m2)

## Decoding states sequence... DONE

Finally, we can see how FPT influences the probability of being in State 1 or State 2:

plotStationary(m2, plotCI=TRUE)

As we’d expect, the model predicts a higher probability of being in State 1 (putative ARS/foraging) when First Passage Time is high. We could do this analysis for any other covariate we had on hand as well.

Behavioral Change Point Analysis (BCPA)

The next method we’re going to take a look at is the behavioral change point analysis (BCPA), which looks for the points in a time series during which there are notable shifts. In our case, we will be applying the method to a movement trajectory to see where an animal may transition between behavioral states, but technically change point analyses can be performed on any time series data (e.g., fluctuating stock values over time or carbon dioxide concentration in the atmosphere over time).

Just as with all other packages, bcpa has its own data format that it prefers, so we will use the bcpa::MakeTrack command to translate a set of 100 coordinates (from our 1326 point albatross path, for the sake of readability in the outputs) into a usable format:

alba_bpca <- bcpa::MakeTrack(albatross10$location_long,albatross10$location_lat,albatross10$timestamp)
plot(alba_bpca)

To obtain the step length and turning angles, use the bcpa::GetVT command, which decomposes the data into single steps and calculates all the statistics:

albatross.VT <- GetVT(alba_bpca)
head(albatross.VT)

##              Z.start              Z.end            S        Phi     Theta
## 2 -89.62097-1.38952i -89.62086-1.38949i 1.055716e-04  0.2254499  2.967982
## 3 -89.62086-1.38949i -89.62098-1.38952i 1.225402e-04 -2.9336341 -3.159084
## 4 -89.62098-1.38952i -89.62092-1.38953i 6.858608e-05 -0.1862436  2.747391
## 5 -89.62092-1.38953i -89.62101-1.38951i 9.487602e-05  2.8585153  3.044759
## 6 -89.62101-1.38951i -89.62095-1.38952i 5.426389e-05 -0.2628578 -3.121373
## 7 -89.62095-1.38952i -89.62098-1.38952i 2.901672e-05  2.9719114  3.234769
##      T.start    T.end     T.mid       dT            V             T.POSIX
## 2 0.05416667 1.532222 0.7931942 1.478055 7.142605e-05 2008-06-23 18:46:19
## 3 1.53222167 3.034999 2.2836106 1.502778 8.154246e-05 2008-06-23 20:15:44
## 4 3.03499944 4.527778 3.7813886 1.492778 4.594525e-05 2008-06-23 21:45:36
## 5 4.52777778 6.036944 5.2823608 1.509166 6.286652e-05 2008-06-23 23:15:40
## 6 6.03694389 7.533333 6.7851386 1.496389 3.626322e-05 2008-06-24 00:45:50
## 7 7.53333333 9.035833 8.2845832 1.502500 1.931230e-05 2008-06-24 02:15:48

The essence of a change point analysis is a sweep across a time series in search of breaks. This sweep can be conducted in a number of ways, but we will focus here on the window sweep, whereby we identify an appropriate windowsize and sensitivity (K) and then the algorithm searches across the time series in search of break points. One can also input a function as the second argument (it can represent any combination of the elements of our albatross.VT dataframe), to serve as a response variable. In this case, we will define a very simple function persistence of movement in a given direction (i.e. ‘persistence velocity’) because we dont really have any a priori conception of what exactly causes change points in this path.

albatross.ws <- WindowSweep(albatross.VT, "V*cos(Theta)", windowsize=100, progress=FALSE, K=2)

The object that is returned by this function (which takes a little while to run, hence our reduction of the dataset to a smaller length) is a ws data frame whose final column indicates proposed break points should be and the parameter values associated with before and after those break point.

head(albatross.ws$ws)

##   Model       LL       bic           mu1           s1      rho1        mu2
## 1     7 720.7501 -1409.194 -2.458589e-05 3.673670e-05 1.0837768 0.11439631
## 2     7 720.7322 -1409.159 -2.449844e-05 3.663456e-05 0.9758211 0.11722744
## 3     7 712.8886 -1393.471 -2.377652e-05 3.629145e-05 0.9474522 0.11278957
## 4     7 705.5419 -1378.778 -2.353749e-05 3.646370e-05 0.9346491 0.10748713
## 5     7 697.9346 -1363.563 -2.303054e-05 3.642509e-05 0.9205235 0.10344379
## 6     7 690.3318 -1348.358 -2.285653e-05 3.663489e-05 0.9218635 0.09985845
##          s2     rho2 Break.bb.time
## 1 0.1337935 7.588469      119.2833
## 2 0.1318655 7.255237      120.7824
## 3 0.1303602 7.361823      120.7824
## 4 0.1298768 7.631846      120.7824
## 5 0.1285572 7.821147      120.7824
## 6 0.1271208 7.999981      120.7824

We can take a look at these suggested breakpoints by looking at the smoothed plot (i.e., the summary in which all the windows are averaged to obtain the “smooth” model). In this plot, the vertical lines represent the significant change points, the width of the lines is proportional to the number of time that change point was selected.

plot(albatross.ws, type="smooth", mu.where="topright")

That doesn’t offer the clearest picture. We can see that there are many change points that have some support. We could, however, add a threshold parameter, which indicates how many of the windows that were swept over the data must have selected a particular changepoint for it to be considered significant. Here, we will use 20 and see what it looks like:

plot(albatross.ws, type="smooth", threshold=20, mu.where="topright")

This reduces our number of change points down to a more reasonable 8, and all of them appear to signify reasonable shifts in our response variable (which combines velocity and angle).

An alternative way to search for change points is to use the ‘flat’ rather than ‘smooth’ method. This analysis first selects changepoints that it deems significant by clustering neighboring change points, and then estimates a homogeneous behavior that occurs between those changepoints.

plot(albatross.ws, type="flat", mu.where="topright")

Once again, if we don’t set an equivalent to the threshold parameter (in the case of the ‘flat’ approach, its called clusterwidth), we get quite a few change points. If we set this parameter to 20, we get the following:

plot(albatross.ws, type="flat", clusterwidth=20, mu.where="topright")

This fairly conservative approach results in 19 significant change points in our time series. A summary of these change points can be obtained using the bcpa::ChangePointSummary command. Here ‘mu.hat’, ‘s.hat’, and ‘rho.hat’ are the estimated mean, standard deviation, and autocorrelation of our response variable (velocity persistence) in each phase.

summary <- ChangePointSummary(albatross.ws, clusterwidth=20)
head(summary$phases)

##          t.cut        mu.hat        s.hat    rho.hat          t0       t1
## 1 (-0.207,120] -2.458589e-05 3.673670e-05 1.08377685  -0.2068058 120.0707
## 2    (120,237]  9.173512e-02 1.286393e-01 4.86922976 120.0706708 236.8816
## 3    (237,392]  5.350859e-02 9.301482e-02 1.89849136 236.8815555 392.3500
## 4    (392,496]  8.122618e-03 4.498000e-02 0.03731658 392.3500292 495.9869
## 5    (496,560]  1.855522e-01 1.919126e-01 2.85542760 495.9868900 559.6516
## 6    (560,692] -1.945615e-05 7.206794e-05 1.08645573 559.6516260 692.4168
##    interval
## 1 120.27748
## 2 116.81088
## 3 155.46847
## 4 103.63686
## 5  63.66474
## 6 132.76521

This summmary suggests eight phases. We can also visualize the path itself with the associated change points using the bcpa::PathPlot command or the bcpa::PhasePlot command:

par(mfrow=c(1,2))
PathPlot(alba_bpca, albatross.ws, type="flat", clusterwidth = 20, main="Flat BCPA", plotlegend=TRUE)
PathPlot(alba_bpca, albatross.ws, type="smooth", main="Smooth BCPA", plotlegend=TRUE)

We see similar patterns as in our previous methods, with periods of high direction persistence during transiting between breeding/foraging areas.